function Nu = initialize_random_environment

    % This script is meant to exactly replicate the values of the moments
    % in the first version of the paper, as they are affected by the bootstrap. 
    % The structure Nu actually used in the analysis is re-initialized later on.

    rng(1005);
    Nu.NumberOfDistricts        = 512-0;
    Rho.Deltat                  = 1/10;                                                            % Size of time step, expressed in fractions of a month (1/4 = one week, 1/30 = one day, etc)
    Nu.NExp                     = 5e1+1;                                                           % Number of sample paths using which expectations are computed 
    Nu.HorizonExp               = 36;                                                              % Horizons of truncation for expectations (in months)
    Nu.TExp                     = floor(Nu.HorizonExp/Rho.Deltat);
    Nu.epsilonExp               = randn(Nu.TExp,Nu.NExp);                                          % Shocks to generate the sample paths for computing expectations
    Nu.epsilonExp(1,:)          = 0;
    Nu.SSim                     = 50;                                                              % Number of simulations used for estimation
    Nu.NSim                     = Nu.NumberOfDistricts;                                                           % Number of districts simulated 
    Nu.HorizonMaxSim            = 20*12;                                                           % Total ength of individual simulation (in months)
    Nu.TMaxSim                  = floor(Nu.HorizonMaxSim/Rho.Deltat);
    Nu.HorizonSim               = 1*14;                                                            % Observation window for each district (in months)
    Nu.TSim                     = floor(Nu.HorizonSim/Rho.Deltat);
    Nu.HorizonShock             = 20*12-8;                                                         % Date at which the aggregate shock occurs
    Nu.DateShock                = floor(Nu.HorizonShock/Rho.Deltat)+2;
    Nu.SSimDiagnostics          = 1+Nu.SSim;                                                       % Number of simulations used in diagnostics; must be larger than Nu.SSim;
    for s=1:Nu.SSimDiagnostics
        epsilonSim{s}           = randn(Nu.TMaxSim,Nu.NSim);                                       % Shocks for simulations
        epsilonSim{s}(1,:)      = 0;
    end
    for s=1:Nu.SSimDiagnostics
        etaSim{s}               = randn(1,Nu.NSim);                                                % Shocks for simulations
    end
    Nu.ainit                    = binornd(1,0.5*ones(1,Nu.NSim));                                  % Initial adoption decisio5
    Nu.Xinit                    = rand(1,Nu.NSim);                                                 % Initial participation rate   
    Nu.B                        = 1e2;
    Nu.Nregions                 = Nu.NumberOfDistricts;
    Nu.Nperiods                 = 13;
    for b=1:Nu.B
        Nu.indices1(:,b)        = randsample(Nu.Nregions,Nu.Nregions,true);
        Nu.indices3(:,b)        = randsample(Nu.Nperiods-5,Nu.Nperiods-5,true);
        Nu.indices4(:,b)        = randsample(Nu.Nregions,Nu.Nregions,true);
    end

end
